To quantitatively examine the efficacy of vegetation restoration in drylands globally.
Study-level viz to document patterns in exclusions primarily and the relatie frequenices, at the study level, of major categories of evidence.
#study data####
library(tidyverse)
studies <- read_csv("data/studies.csv")
studies
## # A tibble: 278 x 18
## ID title technique data region exclude rationale observations
## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 152 Shor~ seeding,~ expe~ Africa no <NA> <NA>
## 2 180 Rest~ chemical~ App.~ Africa no <NA> <NA>
## 3 229 Infl~ soil see~ fiel~ Africa no <NA> <NA>
## 4 230 Acti~ planting fiel~ Africa no <NA> <NA>
## 5 255 The ~ grazing ~ fiel~ Africa no <NA> <NA>
## 6 262 Reve~ seeding,~ eper~ Africa no <NA> <NA>
## 7 263 The ~ phytogen~ fiel~ Africa no <NA> <NA>
## 8 264 Eval~ seeding,~ fiel~ Africa no <NA> <NA>
## 9 271 Patc~ natural ~ fiel~ Africa no <NA> <NA>
## 10 4 Fact~ natural ~ App.~ Africa no <NA> <NA>
## # ... with 268 more rows, and 10 more variables: disturbance <chr>,
## # system <chr>, goal <chr>, intervention <chr>, paradigm <chr>,
## # grazing <chr>, hypothesis <chr>, soil <chr>, benchmark <chr>,
## # notes <chr>
#quick look at rationale needed
exclusions <- studies %>%
filter(exclude == "yes")
#quick look at studies with paradigms
evidence <- studies %>%
filter(exclude == "no")
#library(skimr)
#skim(evidence)
#study-level viz#####
#exclusions
#ggplot(exclusions, aes(rationale, fill = region)) +
# geom_bar() +
# coord_flip() +
# labs(x = "rational for exclusion", y = "frequency") +
# scale_fill_brewer(palette = "Paired")
ggplot(evidence, aes(disturbance, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
#ggplot(evidence, aes(region, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
# coord_flip() +
# scale_fill_brewer(palette = "Paired") +
# labs(y = "frequency")
#ggplot(evidence, aes(data, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
#coord_flip() +
#scale_fill_brewer(palette = "Paired") +
#labs(y = "frequency")
#ggplot(evidence, aes(system, fill = paradigm)) +
# geom_bar(na.rm = TRUE) +
# coord_flip() +
# scale_fill_brewer(palette = "Paired") +
# labs(y = "frequency")
#ggplot(evidence, aes(goal, fill = paradigm)) +
#geom_bar(na.rm = TRUE) +
#coord_flip() +
#scale_fill_brewer(palette = "Paired") +
#labs(x = "outcome", y = "frequency")
#step 1 models####
#paradigm
derived.evidence <- evidence %>%
group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())
#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## region 6 0.301367 160 9.5682 0.9995
## paradigm:region 6 0.213627 154 9.3546 0.9998
#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## goal 6 0.240941 160 9.6287 0.9997
## paradigm:goal 4 0.301480 156 9.3272 0.9897
#even split between active and passive evidence by all key categories
A summary of sort process using PRISMA.
library(PRISMAstatement)
prisma(found = 1504,
found_other = 5,
no_dupes = 1039,
screened = 1039,
screen_exclusions = 861,
full_text = 178,
full_text_exclusions = 101,
qualitative = 77,
quantitative = 66,
width = 800, height = 800)
Check data and calculate necessary measures.
#all data
#data <- read_csv("data/data.csv")
#data <- data %>%
#mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))
#other effect size estimates
#library(compute.es)
#data <- data %>%
#mutate(d=mes(mean.t, mean.c, sd.t, sd.c, n.t, n.c, level = 95, #cer = 0.2, dig = 2, , id = ID, data = data))
#check metafor
#data from ag and grazing studies that examined restoration in drylands
data <- read_csv("data/data.csv")
mydata<- data %>%
filter(disturbance %in% c("agriculture","grazing")) %>%
filter(!notes %in% "couldnt extract data")
#write.csv(mydata, file="mydata.csv")
mydata <- mydata %>%
mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/(n.t*mean.t^2)) + (sd.c^2/(n.c*mean.c^2)))) #use only lrr now
#mydata<- mydata %>%
# mutate(aridity.range = cut(aridity.index, breaks=c(0,5,10,20,25,69.92), labels=c("hyperarid","arid", #"semiarid", "moderately arid", "slightly humid"))) #categorization of aridity.index values, according to de Martone
Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.
#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
#map + geom_point(data=paperdata, aes(x=long, y=lat)) +
#labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally
#add in levels and color code points on map####
#map of 178 articles
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) +
scale_color_brewer(palette = "Paired") +
labs(x = "longitude", y = "latitude", color = "")
#map of selected articles, agriculture and grazing disturbances
map + geom_point(data=mydata, aes(x=long, y=lat, color = paradigm)) +
scale_color_brewer(palette = "Paired") +
labs(x = "longitude", y = "latitude", color = "")
#aggregation####
se <- function(x){
sd(x)/sqrt(length(x))
}
data.simple <- mydata %>%
group_by(study.ID, paradigm, technique, measure.success) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
main.data <- mydata %>%
group_by(study.ID, paradigm, intervention, outcome) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
#EDA data####
simple.data <- mydata %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)
parad.data <- mydata %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)
tech.data <- mydata %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)
success.data <- mydata %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)
#active
active <- mydata %>%
filter(paradigm == "active")
#viz for aggregation####
disturbance.data <- mydata %>%
group_by(study.ID,disturbance, paradigm) %>%
count()
disturbance.data2 <- disturbance.data %>%
group_by(disturbance,paradigm) %>%
count()
ggplot(na.omit(disturbance.data2), aes(disturbance,nn, fill=paradigm))+
geom_bar(stat = "identity") +
coord_flip(ylim=0:44) +
scale_fill_brewer(palette = "Blues")
intervention.data <- active %>%
group_by(study.ID,intervention, paradigm) %>%
count()
intervention.data2 <- intervention.data %>%
group_by(intervention,paradigm) %>%
count()
#ggplot(na.omit(intervention.data2), aes(intervention,nn, fill=paradigm)) +
#geom_bar(stat = "identity") +
#coord_flip(ylim=0:44) +
#scale_fill_brewer(palette = "Blues")
technique.data <- mydata %>%
group_by(study.ID,technique, paradigm) %>%
count()
technique.data2 <- technique.data %>%
group_by(technique,paradigm) %>%
count()
technique.data3 <- technique.data2 %>%
group_by(technique) %>%
count()
aridity<- mydata %>%
group_by(continent) %>% summarize(mean.aridity=mean(aridity.index))
aridity
## # A tibble: 6 x 2
## continent mean.aridity
## <chr> <dbl>
## 1 Africa 9.04
## 2 Asia 20.7
## 3 Europe 26.8
## 4 North America 19.6
## 5 Oceania 16.4
## 6 South America 22.9
#table(mydata$aridity.index)
ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
ggplot(main.data, aes(outcome, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
Exploratory data analyses to understand data and QA/QC using Rii.
Meta and conventional statistical models to explore relative efficacy.
Step 5. Synthesis stats
#p-value meta
#library(metap)
#all data for metas but cleaned for na's
#all data for meta cleaned
mdata.all <- mydata %>%
filter(!is.na(lrr)) %>%
filter(!is.na(var.es)) %>%
filter(!is.na(n.t)) %>%
filter(!is.na(p)) %>%
filter(!is.na(intervention)) %>%
filter(is.finite(lrr)) %>%
filter(!is.na(exp.length)) %>%
filter(!is.na(MAP)) %>%
filter(!is.na(aridity.index)) #%>%
#filter(!is.infinite(var.es))
summary(mdata.all$var.es)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00076 0.15430 0.01441 39.06266
plot(density(mdata.all$var.es))
top_n(mdata.all, -10, var.es)
## # A tibble: 10 x 50
## study.ID ID publication.year data.year exp.year exp.length
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 111 111. 2016 2003 2003 132
## 2 111 111. 2016 2003 2003 132
## 3 111 111. 2016 2003 2003 132
## 4 111 111. 2016 2003 2003 132
## 5 111 111. 2016 2003 2003 132
## 6 111 111. 2016 2003 2003 132
## 7 111 111. 2016 2003 2003 132
## 8 111 111. 2016 2003 2003 132
## 9 111 111. 2016 2003 2003 132
## 10 111 111. 2016 2003 2003 132
## # ... with 44 more variables: disturbance <chr>, focus <chr>,
## # technique <chr>, intervention <chr>, paradigm <chr>, hypothesis <chr>,
## # pathway <chr>, plant.species <chr>, target.plant <chr>,
## # measure.success <chr>, outcome <chr>, measured.factor <chr>,
## # treatment <chr>, control <chr>, unit <chr>, Nsites <dbl>, n.t <dbl>,
## # n.c <dbl>, ntotalsamples <dbl>, mean.t <dbl>, mean.c <dbl>,
## # sd.t <dbl>, sd.c <dbl>, se.t <dbl>, se.c <dbl>, p <dbl>, df <dbl>,
## # measure.dispersion <chr>, lat <dbl>, long <dbl>, continent <chr>,
## # country <chr>, ecosystem <chr>, elevation <dbl>, MAP <dbl>, MAT <dbl>,
## # aridity.index <dbl>, `potential.evaporation (mm)` <dbl>,
## # grazing <chr>, soil <chr>, notes <chr>, lrr <dbl>, rii <dbl>,
## # var.es <dbl>
top_n(mdata.all, 5, var.es)
## # A tibble: 5 x 50
## study.ID ID publication.year data.year exp.year exp.length disturbance
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 11 11.7 2015 2007 1956 600 agriculture
## 2 11 11.7 2015 2007 1956 600 agriculture
## 3 11 11.7 2015 2007 1956 600 agriculture
## 4 11 11.7 2015 2007 1956 600 agriculture
## 5 147 147. 2012 2005 2005 36 agriculture
## # ... with 43 more variables: focus <chr>, technique <chr>,
## # intervention <chr>, paradigm <chr>, hypothesis <chr>, pathway <chr>,
## # plant.species <chr>, target.plant <chr>, measure.success <chr>,
## # outcome <chr>, measured.factor <chr>, treatment <chr>, control <chr>,
## # unit <chr>, Nsites <dbl>, n.t <dbl>, n.c <dbl>, ntotalsamples <dbl>,
## # mean.t <dbl>, mean.c <dbl>, sd.t <dbl>, sd.c <dbl>, se.t <dbl>,
## # se.c <dbl>, p <dbl>, df <dbl>, measure.dispersion <chr>, lat <dbl>,
## # long <dbl>, continent <chr>, country <chr>, ecosystem <chr>,
## # elevation <dbl>, MAP <dbl>, MAT <dbl>, aridity.index <dbl>,
## # `potential.evaporation (mm)` <dbl>, grazing <chr>, soil <chr>,
## # notes <chr>, lrr <dbl>, rii <dbl>, var.es <dbl>
library(metafor)
#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)
#metadata <- metadata %>%
# filter(!is.na(LRR)) %>%
# filter(!is.na(LRR_var)) %>%
# filter(!is.na(n.t)) %>%
# filter(!is.na(p)) %>%
# filter(!is.na(intervention)) %>%
# filter(is.finite(lrr)) %>%
# filter(!is.na(exp.length)) %>%
# filter(!is.na(aridity.index))
#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)
max(mdata.all$var.es)/min(mdata.all$var.es)
## [1] 4.800933e+14
#active only data for meta cleaned up
mdata <- mydata %>%
filter(paradigm == "active") %>%
filter(!is.na(lrr)) %>%
filter(!is.na(var.es)) %>%
filter(!is.na(n.t)) %>%
filter(!is.na(p)) %>%
filter(!is.na(intervention)) %>%
filter(is.finite(lrr)) %>%
filter(!is.na(exp.length))
#aggregated data for metas var estimated with central tendency
#note - could also bootstrap mean variance here instead of arithematic mean
simple.mdata <- mdata %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdata.2 <- mdata %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))
simple.mdata.var <- mdata %>%
group_by(intervention) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
simple.mdata2.var <- mdata %>%
group_by(intervention, outcome) %>%
summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))
#metas with p-values####
#schweder(mdata$p)
#sumz(p, data = mdata)
#mdata %>%
#split(.$intervention) %>%
#purrr::map(~sumz(.$p))
#sumlog(mdata$p)
#metas with meta package on effect size measures####
library(meta)
#all data non-aggregated
#m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, data = mdata)
#summary(m)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)
#t-tests if different from 0
tmu <- function(x){t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)}
mdata.all %>%
split(.$paradigm) %>%
purrr::map(~tmu(.$lrr)) #note this uses arithmetic means not estimated means from random effect models
## $active
##
## One Sample t-test
##
## data: x
## t = 7.6083, df = 1101, p-value = 5.943e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2069479 0.3507825
## sample estimates:
## mean of x
## 0.2788652
##
##
## $passive
##
## One Sample t-test
##
## data: x
## t = -7.5438, df = 357, p-value = 3.824e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.4408271 -0.2585133
## sample estimates:
## mean of x
## -0.3496702
#Propose we only use random effect models because of the diversity of studies, differences in the methods and samples that may introduce heterogeneity
m1 <- metagen(lrr, var.es, studlab = ID, byvar = paradigm, comb.fixed=FALSE, data = mdata.all)
summary(m1) #active is positive and passive is negative so should not mix
## Number of studies combined: k = 1460
##
## 95%-CI z p-value
## Random effects model 0.0766 [0.0654; 0.0879] 13.38 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 184527270028.74 [184527270027.80; 184527270029.68]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 184583923516.49 [184583923515.55; 184583923517.44]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49679407227636228416848646.00 1459 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## paradigm = active 1102 0.2184 [ 0.2055; 0.2314]
## paradigm = passive 358 -0.3413 [-0.3753; -0.3073]
## Q tau^2 I^2
## paradigm = active 49671693556322550580400868.00 0.0450 100.0%
## paradigm = passive 4152232320954931871684.00 0.1047 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 911.23 1 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
#should do a t-test right now of paradigm is diff from 0
metareg(m1,aridity.index)
##
## Mixed-Effects Model (k = 1460; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 33269628721436362276444.00
## R^2 (amount of heterogeneity accounted for): 2.20%
##
## Test for Residual Heterogeneity:
## QE(df = 1458) = 48507118675854218692260048.0000, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 2845.3011, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.7855 0.0144 54.3680 <.0001 0.7571 0.8138 ***
## aridity.index -0.0354 0.0007 -53.3414 <.0001 -0.0367 -0.0341 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#interventions
m2 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, comb.fixed=FALSE, subset = paradigm == "active", data = mdata.all)
summary(m2)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 210054235083.87 [210054235082.78; 210054235084.96]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550580400868.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## intervention = vegetation 779 0.1845 [0.1694; 0.1996]
## intervention = soil 248 0.3128 [0.2990; 0.3265]
## intervention = water addition 75 0.6409 [0.5539; 0.7279]
## Q tau^2 I^2
## intervention = vegetation 48393423300974332080088226.00 0.0443 100.0%
## intervention = soil 97513761686148203152860.00 0.0111 100.0%
## intervention = water addition 31132.18 0.1047 99.8%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 226.58 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
# Meta regression for aridity index (de Martonne 1927) and length of experiments (in months)
metareg(m2, ~ aridity.index)
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 44082670438161540188602.00
## R^2 (amount of heterogeneity accounted for): 2.22%
##
## Test for Residual Heterogeneity:
## QE(df = 1100) = 48490937481977699572460620.0000, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1474.0362, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.7703 0.0158 48.7709 <.0001 0.7393 0.8012 ***
## aridity.index -0.0289 0.0008 -38.3932 <.0001 -0.0304 -0.0274 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
metareg(m2, ~ exp.length )
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value): 0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 44082669939996730852248.00
## R^2 (amount of heterogeneity accounted for): 2.22%
##
## Test for Residual Heterogeneity:
## QE(df = 1100) = 48490936933996402964282222.0000, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 961.4837, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0499 0.0085 5.8663 <.0001 0.0332 0.0666 ***
## exp.length 0.0014 0.0000 31.0078 <.0001 0.0013 0.0015 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m3)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 3420000814.26 [3420000812.64; 3420000815.87]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871684.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## intervention = vegetation 125 0.2654 [ 0.2067; 0.3241]
## intervention = grazing exclusion 29 0.1351 [ 0.0270; 0.2431]
## intervention = soil 204 -0.7583 [-0.8196; -0.6970]
## Q tau^2 I^2
## intervention = vegetation 4152209524073903423668.00 0.1047 100.0%
## intervention = grazing exclusion 238316232.18 0.0881 100.0%
## intervention = soil 14453104123321616.00 0.1990 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 595.91 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
#trying models by continents
m4 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m4)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 210341521299.85 [210341521298.76; 210341521300.95]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550580400868.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## continent = Asia 652 0.2110 [ 0.1947; 0.2273]
## continent = Europe 253 0.3272 [ 0.3136; 0.3408]
## continent = North America 150 0.0531 [ 0.0245; 0.0816]
## continent = South America 19 -0.4895 [-0.6073; -0.3718]
## continent = Africa 16 1.6107 [ 1.3665; 1.8549]
## continent = Oceania 12 -0.0398 [-0.1834; 0.1037]
## Q tau^2 I^2
## continent = Asia 48393423133721090038820246.00 0.0443 100.0%
## continent = Europe 97513761686148203152860.00 0.0111 100.0%
## continent = North America 23018426431807620.00 0.0209 100.0%
## continent = South America 90607.11 0.0643 100.0%
## continent = Africa 254005421361935.34 0.1885 100.0%
## continent = Oceania 1485.24 0.0628 99.3%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 615.00 5 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
m5 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m5)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 3429673787.56 [3429673785.94; 3429673789.18]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871684.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## continent = North America 24 0.1711 [ 0.1162; 0.2260]
## continent = Asia 284 -0.4924 [-0.5448; -0.4400]
## continent = Europe 21 0.0886 [-0.0387; 0.2159]
## continent = South America 4 0.7712 [-0.6724; 2.2148]
## continent = Africa 25 0.3817 [ 0.2547; 0.5087]
## Q tau^2 I^2
## continent = North America 709.76 0.0144 96.8%
## continent = Asia 14453104145902294.00 0.1990 100.0%
## continent = Europe 238286159.01 0.0886 100.0%
## continent = South America 33314949487640844.00 2.1620 100.0%
## continent = Africa 4152172019980625248426.00 0.1047 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 373.50 4 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#use random effects mean and var estimate to plot
#do t-tests here too
#outcomes
m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 210149866430.08 [210149866428.99; 210149866431.17]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550580400868.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## outcome = soil 249 0.2204 [ 0.1558; 0.2849]
## outcome = plants 305 0.5071 [ 0.4936; 0.5206]
## outcome = animals 24 -0.1152 [-0.1155; -0.1148]
## outcome = habitat 524 0.0621 [ 0.0437; 0.0804]
## Q tau^2 I^2
## outcome = soil 35077220764051.67 0.2656 100.0%
## outcome = plants 97513760543782884872406.00 0.0111 100.0%
## outcome = animals 541696.64 <0.0001 100.0%
## outcome = habitat 48393423303339003636244404.00 0.0443 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 8647.81 3 0
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
m <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 3419999089.05 [3419999087.44; 3419999090.67]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871684.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI Q
## outcome = habitat 104 0.1605 [ 0.0964; 0.2246] 4152172019980636258444.00
## outcome = plants 50 0.4438 [ 0.0345; 0.8532] 33314950066906688.00
## outcome = soil 204 -0.7583 [-0.8196; -0.6970] 14453104123321616.00
## tau^2 I^2
## outcome = habitat 0.1047 100.0%
## outcome = plants 2.1620 100.0%
## outcome = soil 0.1990 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 425.72 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
#t-tests for outcomes diff from 0 but you can see using the 95% CI what is different
#check metafor for interactions?? in these big models or are we ok??
#brainstorm on how to explore?? techniques
#save but cut all this.
#m.study <- metagen(lrr, var.es, studlab = study.ID, byvar = intervention, data = mdata)
#summary(m.study)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)
#aggregated data
#m1 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata)
#summary(m1)
#funnel(m1)
#metabias(m1)
#forest(m, sortvar = intervention)
#different var estimate
#m2 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata.var)
#summary(m2)
#funnel(m2)
#metabias(m2)
#forest(m, sortvar = intervention)
#m3 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata.2)
#summary(m3)
#funnel(m)
#radial(m3)
#metabias(m2)
#forest(m, sortvar = intervention)
#m4 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata2.var)
#summary(m4)
#funnel(m)
#radial(m4)
#metabias(m2)
#forest(m, sortvar = intervention)
library(metafor)
# I´ve installed the 'devel' version of metafor (https://wviechtb.github.io/metafor/#installation)
#checking and comparing Lrr and sampling variances values obtained with escalc() function and by functions manually defined (lines 153-155)
#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)
#metadata <- metadata %>%
# filter(!is.na(LRR)) %>%
# filter(!is.na(LRR_var)) %>%
# filter(!is.na(n.t)) %>%
# filter(!is.na(p)) %>%
# filter(!is.na(intervention)) %>%
# filter(is.finite(lrr)) %>%
# filter(!is.na(exp.length)) %>%
# filter(!is.na(aridity.index))
#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)
#write.csv(metadata, "metadata.csv")
#mdata.veg <- mydata %>%
# filter(paradigm == "active") %>%
#filter(intervention=="vegetation") %>% #filtering by intervention
#filter(!is.na(lrr)) %>%
#filter(!is.na(var.es)) %>%
# filter(!is.na(n.t)) %>%
# filter(!is.na(p)) %>%
# filter(!is.na(intervention)) %>%
# filter(is.finite(lrr)) %>%
# filter(!is.na(exp.length))
#mdata.mean<- mdata.all %>%
# group_by(study.ID, aridity.index, exp.length, paradigm, intervention, outcome) %>%
# summarize(mean_lrr= mean(lrr),
# mean_var.es= mean(var.es)) #collapsing data to means
#mdata.mean<- mdata.mean %>%
# filter(!is.infinite(mean_var.es)) #filter out outliers that are infinite
#res1 <- rma(lrr, var.es, mods= aridity.index, data=mdata.all, subset=intervention=="vegetation")
#res2 <- rma(lrr, var.es, data=mdata.all, subset=intervention=="soil") error
mod.1<-rma(lrr,var.es, mods= ~ aridity.index, subset= paradigm == "active", digits=4,data=mdata.all)
mod.1
##
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.8019 (SE = 0.0353)
## tau (square root of estimated tau^2 value): 0.8955
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability): 113092375770.11
## R^2 (amount of heterogeneity accounted for): 8.78%
##
## Test for Residual Heterogeneity:
## QE(df = 1100) = 6499005541922.0801, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 112.0030, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.8574 0.0651 13.1675 <.0001 0.7298 0.9851 ***
## aridity.index -0.0330 0.0031 -10.5831 <.0001 -0.0391 -0.0269 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.2<-rma.mv(yi=lrr, V=var.es, mods= aridity.index, method = "REML", test="t", random = ~ 1 | ID, data=mdata.all, sparse=TRUE)
mod.3 <- rma(yi=lrr, vi=var.es, data = mdata.all)
summary(m3)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 3420000814.26 [3420000812.64; 3420000815.87]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871684.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## intervention = vegetation 125 0.2654 [ 0.2067; 0.3241]
## intervention = grazing exclusion 29 0.1351 [ 0.0270; 0.2431]
## intervention = soil 204 -0.7583 [-0.8196; -0.6970]
## Q tau^2 I^2
## intervention = vegetation 4152209524073903423668.00 0.1047 100.0%
## intervention = grazing exclusion 238316232.18 0.0881 100.0%
## intervention = soil 14453104123321616.00 0.1990 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 595.91 2 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mod.4 <- rma(lrr, var.es, mods= cbind (aridity.index, exp.length), data = mdata.all, subset = paradigm == "active")
summary(m4)
## Number of studies combined: k = 1102
##
## 95%-CI z p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 210341521299.85 [210341521298.76; 210341521300.95]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 49671693556322550580400868.00 1101 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## continent = Asia 652 0.2110 [ 0.1947; 0.2273]
## continent = Europe 253 0.3272 [ 0.3136; 0.3408]
## continent = North America 150 0.0531 [ 0.0245; 0.0816]
## continent = South America 19 -0.4895 [-0.6073; -0.3718]
## continent = Africa 16 1.6107 [ 1.3665; 1.8549]
## continent = Oceania 12 -0.0398 [-0.1834; 0.1037]
## Q tau^2 I^2
## continent = Asia 48393423133721090038820246.00 0.0443 100.0%
## continent = Europe 97513761686148203152860.00 0.0111 100.0%
## continent = North America 23018426431807620.00 0.0209 100.0%
## continent = South America 90607.11 0.0643 100.0%
## continent = Africa 254005421361935.34 0.1885 100.0%
## continent = Oceania 1485.24 0.0628 99.3%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 615.00 5 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mod.5 <- rma(lrr, var.es, mods= ~intervention, data = mdata.all, subset = paradigm == "active")
summary(m5)
## Number of studies combined: k = 358
##
## 95%-CI z p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
##
## Quantifying residual heterogeneity:
## H = 3429673787.56 [3429673785.94; 3429673789.18]; I^2 = 100.0% [100.0%; 100.0%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 4152232320954931871684.00 357 0
##
## Results for subgroups (random effects model):
## k 95%-CI
## continent = North America 24 0.1711 [ 0.1162; 0.2260]
## continent = Asia 284 -0.4924 [-0.5448; -0.4400]
## continent = Europe 21 0.0886 [-0.0387; 0.2159]
## continent = South America 4 0.7712 [-0.6724; 2.2148]
## continent = Africa 25 0.3817 [ 0.2547; 0.5087]
## Q tau^2 I^2
## continent = North America 709.76 0.0144 96.8%
## continent = Asia 14453104145902294.00 0.1990 100.0%
## continent = Europe 238286159.01 0.0886 100.0%
## continent = South America 33314949487640844.00 2.1620 100.0%
## continent = Africa 4152172019980625248426.00 0.1047 100.0%
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 373.50 4 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#forest(m1, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m2, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m3, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#forest(m4, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))
#t-tests for lrr diff from 0
#mdata %>%
#split(.$intervention) %>%
#purrr::map(~sumz(.$lrr))
#effect sizes plots
#need better ci estimates
#ggplot(simple.mdata, aes(intervention, lrr)) +
# ylim(c(-2,2)) +
# geom_point(position = position_dodge(width = 0.5)) +
# labs(x = "", y = "lrr") +
# coord_flip() +
# geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = position_dodge(width = 0.5)) +
# geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
#ggplot(simple.mdata.2, aes(intervention, lrr, color = outcome)) +
# ylim(c(-2,2)) +
#geom_point(position = position_dodge(width = 0.5)) +
# labs(x = "", y = "lrr", color = "") +
# coord_flip() +
#geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = #position_dodge(width = 0.5)) +
#geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
#ggplot(mdata, aes(lrr, color = intervention)) +
#geom_freqpoly(binwidth = 0.5, size = 2) +
#xlim(c(-10, 10)) +
#labs(x = "lrr", y = "frequency", color = "") +
#geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
#scale_color_brewer()
#ggplot(mdata, aes(lrr, fill = intervention)) +
#geom_dotplot(binwidth = 1) +
# xlim(c(-10, 10)) +
# labs(x = "lrr", y = "frequency", fill = "") +
# geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
# scale_fill_brewer()